# Load Relevant Packages --------------------------------------------------
library(ggplot2)
library(reshape)
library(patchwork)
library(stargazer)
library(numform)
library(marginaleffects)
library(cem)
library(survey)
library(psych)

# Load Data ---------------------------------------------------------------
# Set file path
setwd("")

load("drop_in_ocean_data.rdata")


# Figure 1 ----------------------------------------------------------------
# * No Effect --------------------------------------------------------------
# Dummy Data Frame
d <- data.frame(group = c(rep("Group 1", 7), rep("Group 2", 7)),
                iv = rep(1:7, 2),
                dv = c(rep(7, 7), rep(2, 7)))

# Store plot
pNull <- ggplot(d, aes(x = iv, y = dv, group = group)) +
  geom_line(size = 2, aes(color = group)) +
  scale_color_manual(values = c("black", "grey")) +
  theme_bw() +
  scale_x_continuous(breaks = c(1, 7), labels = c("Fewer", "More")) +
  scale_y_continuous(limits = c(0, 12), breaks = c(0, 12), labels = c("Positive", "Negative")) +
  labs(y = "Carceral State Attitude", x = "Experiences", title = "No Effect") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")



# * Common Effect -----------------------------------------------------------
# Dummy Data Frame
d <- data.frame(group = c(rep("Group 1", 7), rep("Group 2", 7)),
                iv = rep(1:7, 2),
                dv = c(seq(7, 12, length.out = 7),
                       seq(2, 7, length.out = 7)))

# Store plot
pPara <- ggplot(d, aes(x = iv, y = dv, group = group)) +
  geom_line(size = 2, aes(color = group)) +
  scale_color_manual(values = c("black", "grey")) +
  theme_bw() +
  scale_x_continuous(breaks = c(1, 7), labels = c("Fewer", "More")) +
  scale_y_continuous(limits = c(0, 12), breaks = c(0, 12), labels = c("Positive", "Negative")) +
  labs(y = "Carceral State Attitude", x = "Experiences", title = "Common Effect") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")

# * Convergence -------------------------------------------------------------
# Dummy Data Frame
d <- data.frame(group = c(rep("Group 1", 7), rep("Group 2", 7)),
                iv = rep(1:7, 2),
                dv = c(seq(7, 10, length.out = 7),
                       seq(2, 10, length.out = 7)))

# Store plot
pConverg <- ggplot(d, aes(x = iv, y = dv, group = group)) +
  geom_line(size = 2, aes(color = group)) +
  scale_color_manual(values = c("black", "grey")) +
  theme_bw() +
  scale_x_continuous(breaks = c(1, 7), labels = c("Fewer", "More")) +
  scale_y_continuous(limits = c(0, 12), breaks = c(0, 12), labels = c("Positive", "Negative")) +
  labs(y = "Carceral State Attitude", x = "Experiences", title = "Different Effect: Convergence") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")


# * Polarization ------------------------------------------------------------
# Dummy Data Frame
d <- data.frame(group = c(rep("Group 1", 7), rep("Group 2", 7)),
                iv = rep(1:7, 2),
                dv = c(seq(7, 12, length.out = 7),
                       seq(2, 5, length.out = 7)))

# Store plot
pConverg2 <- ggplot(d, aes(x = iv, y = dv, group = group)) +
  geom_line(size = 2, aes(color = group)) +
  scale_color_manual(values = c("black", "grey")) +
  theme_bw() +
  scale_x_continuous(breaks = c(1, 7), labels = c("Fewer", "More")) +
  scale_y_continuous(limits = c(0, 12), breaks = c(0, 12), labels = c("Positive", "Negative")) +
  labs(y = "Carceral State Attitude", x = "Experiences", title = "Different Effect: Polarization") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")

# * Divergence --------------------------------------------------------------
# Dummy Data Frame
d <- data.frame(group = c(rep("Group 1", 7), rep("Group 2", 7)),
                iv = rep(1:7, 2),
                dv = c(seq(5, 10, length.out = 7),
                       seq(2, 1, length.out = 7)))

# Store Plot
pDiverge <- ggplot(d, aes(x = iv, y = dv, group = group)) +
  geom_line(size = 2, aes(color = group)) +
  scale_color_manual(values = c("black", "grey")) +
  theme_bw() +
  scale_x_continuous(breaks = c(1, 7), labels = c("Fewer", "More")) +
  scale_y_continuous(limits = c(0, 12), breaks = c(0, 12), labels = c("Positive", "Negative")) +
  labs(y = "Carceral State Attitude", x = "Experiences", title = "Different Effect: Divergence") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")


# ** Combine Plots and Store -------------------------------------------------
pNull + pPara + pConverg2 + pConverg + #pDiverge +
  plot_layout(guides = "collect") & theme(legend.position = 'bottom') &
  plot_annotation(tag_levels = 'A')
# ggsave("figure1.pdf", height = 6, width = 8)


# Figure 2 ----------------------------------------------------------------
# Store experience frequency distributions
# By Race
direct_prop_bw <- prop.table(svytable(~direct_exp + black, d.bw.stk), 2)
direct_prop_bw <- melt(direct_prop_bw)
direct_prop_bw$black <- factor(direct_prop_bw$black,
                               labels = c("White", "Black"))
direct_prop_bw$var <- "direct"
direct_prop_bw$group <- "Race"
names(direct_prop_bw)[c(1, 2)] <- c("x", "person")

# Gender among Blacks
direct_prop_sex_blk <- prop.table(svytable(~direct_exp + man, subset(d.bw.stk, black == 1)), 2)
direct_prop_sex_blk <- melt(direct_prop_sex_blk)
direct_prop_sex_blk$var <- "direct"
direct_prop_sex_blk$man <- factor(direct_prop_sex_blk$man,
                                  labels = c("Women", "Men"))
direct_prop_sex_blk$group <- "Black"
names(direct_prop_sex_blk)[c(1, 2)] <- c("x", "person")

# Gender among Whites
direct_prop_sex_wht <- prop.table(svytable(~direct_exp + man, subset(d.bw.stk, black == 0)), 2)
direct_prop_sex_wht <- melt(direct_prop_sex_wht)
direct_prop_sex_wht$var <- "direct"
direct_prop_sex_wht$man <- factor(direct_prop_sex_wht$man,
                                  labels = c("Women", "Men"))
direct_prop_sex_wht$group <- "White"
names(direct_prop_sex_wht)[c(1, 2)] <- c("x", "person")

# combine gender
direct_prop_int <- rbind(direct_prop_sex_blk, direct_prop_sex_wht)

#  Information for Plot
direct_prop_int$person2 <- paste(direct_prop_int$group, direct_prop_int$person, sep = "\n")

# Level labels
direct_prop_bw$stack_group <- factor(direct_prop_bw$person,
                                     labels = c("White", "Black"))
direct_prop_int$stack_group <- paste(direct_prop_int$group, direct_prop_int$person)

# Combine proportions
group_df <- rbind(subset(direct_prop_bw, select = c("x", "value", "stack_group")), 
                  subset(direct_prop_int, select = c("x", "value", "stack_group")))
group_df$stack_group <- factor(group_df$stack_group,
                               levels = c("Black", "Black Men", "Black Women", "White",  "White Men", "White Women"))

# data frame with group means for plotting
means <- data.frame(stack_group = unique(levels(group_df$stack_group)),
                    mean = c(svymean(~direct_exp, subset(d.bw.stk, black == 1))[[1]],
                             svymean(~direct_exp, subset(d.bw.stk, black == 1 & man == 1))[[1]],
                             svymean(~direct_exp, subset(d.bw.stk, black == 1 & man == 0))[[1]],
                             svymean(~direct_exp, subset(d.bw.stk, black == 0))[[1]],
                             svymean(~direct_exp, subset(d.bw.stk, black == 0 & man == 1))[[1]],
                             svymean(~direct_exp, subset(d.bw.stk, black == 0 & man == 0))[[1]]))
means$stack_group <- factor(means$stack_group,
                            levels = c("Black", "Black Men", "Black Women", "White",  "White Men", "White Women"))

# plot
ggplot(group_df, aes(x = x, y = value)) +
  geom_vline(data = means, aes(xintercept = mean)) +
  geom_col() +
  facet_wrap(~stack_group, ncol = 3) +
  theme_bw() +
  labs(x = "", y = "Proportion") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 14),
        plot.title = element_text(size = 16),
        legend.title = element_blank(),
        legend.text = element_text(size = 12),
        legend.position = "bottom")
# ggsave("figure2.pdf", width = 8, height = 6)






# DV Scale Internal Consistency -------------------------------------------
# Define variables
vars <- c("felon_benefits", "felon_voting", "cit_osight", "pub_defend", "blm_support",
          "p.race.fair", "p.exces.force", "p.account", "court.fair.all")
# Omega calculations
omega(subset(cjs.df, select = vars))
omega(subset(cjs.df, black == 1, select = vars))
omega(subset(cjs.df, black == 0, select = vars))


# Contact Frequency ----------------------------------------------------------------
# * By Race ---------------------------------------------------------------
# bivariate model
m_cjs_att_bivar <- lm(cjs_att ~ direct_exp*black,
                      data = cjs.df, weights = wts_bw_stk)
summary(m_cjs_att_bivar)

# model with controls
cjs_att_mod_race <- lm(cjs_att ~ direct_exp*black + rr_sc*black + 
                         peer.felony*black + age_sc*black + man*black + 
                         educ_sc*black + inc_sc*black + pid7_dem*black,
                       data = cjs.df, weights = wts_bw_stk)
summary(cjs_att_mod_race)


# * By Gender (Black) -----------------------------------------------------
# bivariate model
m_cjs_att_blk_bivar <- lm(cjs_att ~ direct_exp*man,
                          data = cjs.df, subset = black == 1,  weights = wts_bw_stk)
summary(m_cjs_att_blk_bivar)

# model with controls
cjs_att_mod_gender_blk <- lm(cjs_att ~ direct_exp*man + rr_sc*man +
                               peer.felony*man + age_sc*man + educ_sc*man + inc_sc*man + pid7_dem*man,
                             data = cjs.df, subset = black == 1,  weights = wts_bw_stk)
summary(cjs_att_mod_gender_blk)




# * By Gender (White) -----------------------------------------------------
# bivariate model
m_cjs_att_wht_bivar <- lm(cjs_att ~ direct_exp*man,
                          data = cjs.df, subset = black == 0,  weights = wts_bw_stk)
summary(m_cjs_att_wht_bivar)


# model with controls
cjs_att_mod_gender_wht <- lm(cjs_att ~ direct_exp*man + rr_sc*man +
                               peer.felony*man + age_sc*man + educ_sc*man + inc_sc*man + pid7_dem*man,
                             data = cjs.df, subset = black == 0,  weights = wts_bw_stk)
summary(cjs_att_mod_gender_wht)





# Figure 3 ----------------------------------------------------------------
# * Race Facet ------------------------------------------------------------
# * Figure Data Sets -------------------------------------------------
blk_samp_direct <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 1), na.rm = T)[[1]],
                              man = names(svytable(~man, subset(d.bw.stk, black == 1)))[which.max(svytable(~man, subset(d.bw.stk, black == 1)))],
                              educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 1)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 1)))],
                              inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 1)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 1)))],
                              pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 1)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 1)))],
                              direct_exp = seq(0, 1, length.out = 10),
                              peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 1)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 1)))],
                              rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 1), na.rm = T)[[1]],
                              black = 1,
                              cjs_att = NA)
blk_samp_direct <- as.data.frame(apply(blk_samp_direct, 2, as.numeric))

wht_samp_direct <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 0), na.rm = T)[[1]],
                              man = names(svytable(~man, subset(d.bw.stk, black == 0)))[which.max(svytable(~man, subset(d.bw.stk, black == 0)))],
                              educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 0)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 0)))],
                              inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 0)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 0)))],
                              pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 0)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 0)))],
                              direct_exp = seq(0, 1, length.out = 10),
                              peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 0)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 0)))],
                              rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 0), na.rm = T)[[1]],
                              black = 0,
                              cjs_att = NA)
wht_samp_direct <- as.data.frame(apply(wht_samp_direct, 2, as.numeric))

# * Figure Predictions ----------------------------------------------------------
# bivariate predicted values
blk_direct <- predict(m_cjs_att_bivar, blk_samp_direct, 
                      se.fit = TRUE)
wht_direct <- predict(m_cjs_att_bivar, wht_samp_direct, 
                      se.fit = TRUE)

# save data frame
pred_direct_bi <- data.frame(var = "direct_exp",
                             fit = c(blk_direct$fit, 
                                     wht_direct$fit),
                             se = c(blk_direct$se.fit, 
                                    wht_direct$se.fit),
                             black = c(rep(1, 10), rep(0, 10)),
                             x = c(seq(0, 1, length.out = 10), seq(0, 1, length.out = 10)))

eval_preds_bi <- rbind(pred_direct_bi)
eval_preds_bi$black <- factor(eval_preds_bi$black,
                              labels = c("White", "Black"))


# full model predicted values
blk_direct <- predict(cjs_att_mod_race, blk_samp_direct, 
                      se.fit = TRUE)
wht_direct <- predict(cjs_att_mod_race, wht_samp_direct, 
                      se.fit = TRUE)

# save data frame
pred_direct <- data.frame(var = "direct_exp",
                          fit = c(blk_direct$fit, 
                                  wht_direct$fit),
                          se = c(blk_direct$se.fit, 
                                 wht_direct$se.fit),
                          black = c(rep(1, 10), rep(0, 10)),
                          x = c(seq(0, 1, length.out = 10), seq(0, 1, length.out = 10)))

eval_preds <- rbind(pred_direct)
eval_preds$black <- factor(eval_preds$black,
                           labels = c("White", "Black"))


# * Plot ------------------------------------------------------------------
bw_bivar <- ggplot(eval_preds_bi, aes(x = x, y = fit, group = black)) +
  geom_ribbon(aes(ymin = fit - 1.39*se, ymax = fit + 1.39*se, fill = black),
              alpha = .6) +
  geom_line(size = 1, aes(color = black)) +
  scale_color_manual(values = c("grey", "black")) +
  scale_fill_manual(values = c("grey", "black")) +
  theme_bw() +
  scale_x_continuous(labels = ff_num(zero = 0, digits = 2)) +
  scale_y_continuous(labels = ff_num(zero = 0, digits = 2)) +
  labs(x = "Direct Experiences", y = "Predicted Carceral State Attitude") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")


# * Black Figure Data Sets -------------------------------------------------
man_samp_direct <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 1 & man == 1), na.rm = T)[[1]],
                              black = 1,
                              educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 1 & man == 1)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 1 & man == 1)))],
                              inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 1 & man == 1)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 1 & man == 1)))],
                              pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 1 & man == 1)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 1 & man == 1)))],
                              direct_exp = seq(0, 1, length.out = 10),
                              peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 1 & man == 1)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 1 & man == 1)))],
                              rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 1 & man == 1), na.rm = T)[[1]],
                              man = 1,
                              cjs_att = NA)
man_samp_direct <- as.data.frame(apply(man_samp_direct, 2, as.numeric))

wom_samp_direct <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 1 & man == 0), na.rm = T)[[1]],
                              black = 1,
                              educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 1 & man == 0)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 1 & man == 0)))],
                              inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 1 & man == 0)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 1 & man == 0)))],
                              pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 1 & man == 0)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 1 & man == 0)))],
                              direct_exp = seq(0, 1, length.out = 10),
                              peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 1 & man == 0)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 1 & man == 0)))],
                              rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 1 & man == 0), na.rm = T)[[1]],
                              man = 0,
                              cjs_att = NA)
wom_samp_direct <- as.data.frame(apply(wom_samp_direct, 2, as.numeric))

# * Black Figure Predictions ----------------------------------------------------------
# Bivariate predicted values
man_direct <- predict(m_cjs_att_blk_bivar, man_samp_direct, 
                      se.fit = TRUE)
wom_direct <- predict(m_cjs_att_blk_bivar, wom_samp_direct, 
                      se.fit = TRUE)

# store data frame
pred_direct_bi <- data.frame(var = "direct_exp",
                             fit = c(man_direct$fit, 
                                     wom_direct$fit),
                             se = c(man_direct$se.fit, 
                                    wom_direct$se.fit),
                             man = c(rep(1, 10), rep(0, 10)),
                             x = c(seq(0, 1, length.out = 10), seq(0, 1, length.out = 10)))

eval_preds_blk_bi <- rbind(pred_direct_bi)
eval_preds_blk_bi$man <- factor(eval_preds_blk_bi$man,
                                labels = c("Women", "Men"))

# full model predicted values
man_direct <- predict(cjs_att_mod_gender_blk, man_samp_direct, 
                      se.fit = TRUE)
wom_direct <- predict(cjs_att_mod_gender_blk, wom_samp_direct, 
                      se.fit = TRUE)

# save data frame
pred_direct <- data.frame(var = "direct_exp",
                          fit = c(man_direct$fit, 
                                  wom_direct$fit),
                          se = c(man_direct$se.fit, 
                                 wom_direct$se.fit),
                          man = c(rep(1, 10), rep(0, 10)),
                          x = c(seq(0, 1, length.out = 10), seq(0, 1, length.out = 10)))

eval_preds_blk <- rbind(pred_direct)
eval_preds_blk$man <- factor(eval_preds_blk$man,
                             labels = c("Women", "Men"))


# * White Figure Data Sets -------------------------------------------------
man_samp_direct <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 0 & man == 1), na.rm = T)[[1]],
                              black = 0,
                              educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 0 & man == 1)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 0 & man == 1)))],
                              inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 0 & man == 1)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 0 & man == 1)))],
                              pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 0 & man == 1)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 0 & man == 1)))],
                              direct_exp = seq(0, 1, length.out = 10),
                              peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 0 & man == 1)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 0 & man == 1)))],
                              rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 0 & man == 1), na.rm = T)[[1]],
                              man = 1,
                              cjs_att = NA)
man_samp_direct <- as.data.frame(apply(man_samp_direct, 2, as.numeric))

wom_samp_direct <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 0 & man == 0), na.rm = T)[[1]],
                              black = 0,
                              educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 0 & man == 0)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 0 & man == 0)))],
                              inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 0 & man == 0)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 0 & man == 0)))],
                              pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 0 & man == 0)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 0 & man == 0)))],
                              direct_exp = seq(0, 1, length.out = 10),
                              peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 0 & man == 0)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 0 & man == 0)))],
                              rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 0 & man == 0), na.rm = T)[[1]],
                              man = 0,
                              cjs_att = NA)
wom_samp_direct <- as.data.frame(apply(wom_samp_direct, 2, as.numeric))

# * White Figure Predictions ----------------------------------------------------------
# control predictions
man_direct <- predict(cjs_att_mod_gender_wht, man_samp_direct, 
                      se.fit = TRUE)
wom_direct <- predict(cjs_att_mod_gender_wht, wom_samp_direct, 
                      se.fit = TRUE)

# store in data frame
pred_direct <- data.frame(var = "direct_exp",
                          fit = c(man_direct$fit, 
                                  wom_direct$fit),
                          se = c(man_direct$se.fit, 
                                 wom_direct$se.fit),
                          man = c(rep(1, 10), rep(0, 10)),
                          x = c(seq(0, 1, length.out = 10), seq(0, 1, length.out = 10)))

eval_preds_wht <- rbind(pred_direct)
eval_preds_wht$man <- factor(eval_preds_wht$man,
                             labels = c("Women", "Men"))

# bivariate predictions
man_direct <- predict(m_cjs_att_wht_bivar, man_samp_direct, 
                      se.fit = TRUE)
wom_direct <- predict(m_cjs_att_wht_bivar, wom_samp_direct, 
                      se.fit = TRUE)

# store in data frame
pred_direct_bi <- data.frame(var = "direct_exp",
                             fit = c(man_direct$fit, 
                                     wom_direct$fit),
                             se = c(man_direct$se.fit, 
                                    wom_direct$se.fit),
                             man = c(rep(1, 10), rep(0, 10)),
                             x = c(seq(0, 1, length.out = 10), seq(0, 1, length.out = 10)))

eval_preds_wht_bi <- rbind(pred_direct_bi)
eval_preds_wht_bi$man <- factor(eval_preds_wht_bi$man,
                                labels = c("Women", "Men"))



# * Gender Facet ---------------------------------------------------
# combine predictions for gender
eval_preds_wht_bi$race <- "White"
eval_preds_blk_bi$race <- "Black"

preds_comb_bi <- rbind(eval_preds_wht_bi, eval_preds_blk_bi)


# Figure
gend_bivar <- ggplot(preds_comb_bi, aes(x = x, y = fit, group = man)) +
  geom_ribbon(aes(ymin = fit - 1.39*se, ymax = fit + 1.39*se, fill = man),
              alpha = .6) +
  geom_line(size = 1, aes(color = man)) +
  scale_color_manual(values = c("grey", "black")) +
  scale_fill_manual(values = c("grey", "black")) +
  scale_x_continuous(labels = ff_num(zero = 0, digits = 2)) +
  scale_y_continuous(labels = ff_num(zero = 0, digits = 2)) +
  facet_grid(~race) +
  theme_bw() +
  labs(x = "Direct Experiences", y = "Predicted Carceral State Attitude") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")


# * Combine and Save --------------------------------------------------------

bw_bivar/gend_bivar + plot_annotation(tag_levels = 'A')
# ggsave("figure3.pdf", width = 6, height = 8)


# Figure 4 ----------------------------------------------------------------
# Figure
ggplot(eval_preds, aes(x = x, y = fit, group = black)) +
  geom_ribbon(aes(ymin = fit - 1.39*se, ymax = fit + 1.39*se, fill = black),
              alpha = .6) +
  geom_line(size = 1, aes(color = black)) +
  scale_color_manual(values = c("grey", "black")) +
  scale_fill_manual(values = c("grey", "black")) +
  theme_bw() +
  scale_x_continuous(labels = ff_num(zero = 0, digits = 2)) +
  scale_y_continuous(labels = ff_num(zero = 0, digits = 2)) +
  labs(x = "Direct Experiences", y = "Predicted Carceral State Attitude") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")
# ggsave("figure4.pdf", width = 6, height = 4)



# Figure 5 ----------------------------------------------------------------
# * Plot --------------------------------------------------------------------
# combine predictions for plot
eval_preds_wht$race <- "White"
eval_preds_blk$race <- "Black"

preds_comb <- rbind(eval_preds_wht, eval_preds_blk)


# Figure
ggplot(preds_comb, aes(x = x, y = fit, group = man)) +
  geom_ribbon(aes(ymin = fit - 1.39*se, ymax = fit + 1.39*se, fill = man),
              alpha = .6) +
  geom_line(size = 1, aes(color = man)) +
  scale_color_manual(values = c("grey", "black")) +
  scale_fill_manual(values = c("grey", "black")) +
  scale_x_continuous(labels = ff_num(zero = 0, digits = 2)) +
  scale_y_continuous(labels = ff_num(zero = 0, digits = 2)) +
  facet_grid(~race) +
  theme_bw() +
  labs(x = "Direct Experiences", y = "Predicted Carceral State Attitude") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")
# ggsave("figure5.pdf", width = 6, height = 4)

# Figure 6 ----------------------------------------------------------------
# Store experience quality distributions
# By Race
direct_Qual_prop_bw <- prop.table(svytable(~pol_quality + black, d.bw.stk), 2)
direct_Qual_prop_bw <- melt(direct_Qual_prop_bw)
direct_Qual_prop_bw$black <- factor(direct_Qual_prop_bw$black,
                                    labels = c("White", "Black"))
direct_Qual_prop_bw$var <- "quality"
direct_Qual_prop_bw$group <- "Race"
names(direct_Qual_prop_bw)[c(1, 2)] <- c("x", "person")

# by gender among Blacks
direct_Qual_prop_sex_blk <- prop.table(svytable(~pol_quality + man, subset(d.bw.stk, black == 1)), 2)
direct_Qual_prop_sex_blk <- melt(direct_Qual_prop_sex_blk)
direct_Qual_prop_sex_blk$var <- "quality"
direct_Qual_prop_sex_blk$man <- factor(direct_Qual_prop_sex_blk$man,
                                       labels = c("Women", "Men"))
direct_Qual_prop_sex_blk$group <- "Black"
names(direct_Qual_prop_sex_blk)[c(1, 2)] <- c("x", "person")

# by gender among Whites
direct_Qual_prop_sex_wht <- prop.table(svytable(~pol_quality + man, subset(d.bw.stk, black == 0)), 2)
direct_Qual_prop_sex_wht <- melt(direct_Qual_prop_sex_wht)
direct_Qual_prop_sex_wht$var <- "quality"
direct_Qual_prop_sex_wht$man <- factor(direct_Qual_prop_sex_wht$man,
                                       labels = c("Women", "Men"))
direct_Qual_prop_sex_wht$group <- "White"
names(direct_Qual_prop_sex_wht)[c(1, 2)] <- c("x", "person")

# combine gender
direct_Qual_prop_int <- rbind(direct_Qual_prop_sex_blk, direct_Qual_prop_sex_wht)

#  Information for Plot
direct_Qual_prop_int$person2 <- paste(direct_Qual_prop_int$group, direct_Qual_prop_int$person, sep = "\n")

direct_Qual_prop_bw$x <- factor(direct_Qual_prop_bw$x,
                                labels = c("Safe and Fair", "Unsafe or Unfair", "Unsafe and Unfair"))
direct_Qual_prop_int$x <- factor(direct_Qual_prop_int$x,
                                 labels = c("Safe and Fair", "Unsafe or Unfair", "Unsafe and Unfair"))

# labels
direct_Qual_prop_bw$stack_group <- factor(direct_Qual_prop_bw$person,
                                          labels = c("White", "Black"))
direct_Qual_prop_int$stack_group <- paste(direct_Qual_prop_int$group, direct_Qual_prop_int$person)

# Combine to plot
group_df <- rbind(subset(direct_Qual_prop_bw, select = c("x", "value", "stack_group")), 
                  subset(direct_Qual_prop_int, select = c("x", "value", "stack_group")))

# Relabel for plot
group_df$stack_group <- factor(group_df$stack_group,
                               levels = c("Black", "Black Men", "Black Women", "White",  "White Men", "White Women"))
group_df$x2 <- NA
group_df$x2[which(group_df$x == "Safe and Fair")] <- "Positive"
group_df$x2[which(group_df$x == "Unsafe or Unfair")] <- "Mixed"
group_df$x2[which(group_df$x == "Unsafe and Unfair")] <- "Negative"
group_df$x2 <- factor(group_df$x2,
                      levels = c("Positive", "Mixed", "Negative"))

# Plot
ggplot(group_df, aes(x = x2, y = value)) +
  geom_col() +
  facet_wrap(~stack_group, ncol = 3) +
  theme_bw() +
  labs(x = "", y = "Proportion") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title.y = element_text(size = 14),
        plot.title = element_text(size = 16),
        legend.title = element_blank(),
        legend.text = element_text(size = 12),
        legend.position = "bottom")
# ggsave("figure6.pdf", width = 8, height = 6)



# Contact Quality ----------------------------------------------------------------
# * Create Quality IV -------------------------------------------------------
cjs.df$pol_quality <- factor(cjs.df$pol_quality,
                             labels = c("Positive", "Mixed", "Negative"))


# * By Race -----------------------------------------------------------------
# ** Specify Matching Weights ------------------------------------------------
# define variables
demogs <- c("pid3", "rr_sc", "peer.felony", "age", "man", "educ3", "inc")

# subset data
cem_dat <- na.omit(subset(cjs.df, 
                          select = c(demogs, "pol.stop","pol_quality", "cjs_att", "black", "wts_bw_stk",
                                     "pid7_dem", "age_sc", "educ_sc", "inc_sc")))

# define breaks for continuous items
cuts <- list(age = 4, inc = c(-.05, 2, 6, 11), peer.felony = c(-.05, .33, 1),
             rr_sc = 4) 

# create matched object
mat <- cem(treatment = "pol_quality", 
           data = cem_dat, 
           drop = c("cjs_att", "black", "wts_bw_stk",
                    "pid7_dem", "age_sc", "educ_sc", "inc_sc"),
           cutpoints = cuts)
mat


# ** Estimate Model ----------------------------------------------------------
m_race <- lm(cjs_att ~ pol_quality*black + pol.stop*black + rr_sc*black + peer.felony*black + age_sc*black + 
               man*black + educ_sc*black + inc_sc*black + pid7_dem*black, 
             cem_dat, weights = mat$w)
summary(m_race)


# * By Gender within Race -----------------------------------------------------------------
# ** Specify Matching Weights ------------------------------------------------
# specify variables
demogs <- c("pid3", "rr_sc", "peer.felony", "man", "age", "educ3", "inc")

# Black respondents subset
cem_dat_black <- na.omit(subset(cjs.df,
                                black == 1,
                                select = c(demogs, "pol.stop","pol_quality", "cjs_att", "wts_bw_stk",
                                           "pid7_dem", "age_sc", "educ_sc", "inc_sc")))

# define breaks for continuous variables
cuts_blk <- list(age = c(17, 28, 40, 90), inc = c(-.05, 2, 6, 11), peer.felony = c(-.05, .33, 1), 
                 rr_sc = 4) 

# create matched object
mat_black <- cem(treatment = "pol_quality", 
                 data = cem_dat_black, 
                 drop = c("cjs_att", "man", "wts_bw_stk",
                          "pid7_dem", "age_sc", "educ_sc", "inc_sc"),
                 cutpoints = cuts_blk
)
mat_black

# White respondents subset
cem_dat_white <- na.omit(subset(cjs.df,
                                black == 0,
                                select = c(demogs, "pol.stop","pol_quality", "cjs_att", "wts_bw_stk",
                                           "pid7_dem", "age_sc", "educ_sc", "inc_sc")))

# define breaks for continuous variables
cuts_wht <- list(age = c(17, 28, 38, 56, 90), inc = c(-.05, 2, 5, 11), peer.felony = 2, 
                 rr_sc = 4) 

# create matched object
mat_white <- cem(treatment = "pol_quality", 
                 data = cem_dat_white, 
                 drop = c("cjs_att", "man", "wts_bw_stk",
                          "pid7_dem", "age_sc", "educ_sc", "inc_sc"),
                 cutpoints = cuts_wht
)
mat_white


# ** Estimate Model ----------------------------------------------------------
# Black respondents
m_gend_blk <- lm(cjs_att ~ pol_quality*man + pol.stop*man + rr_sc*man + peer.felony*man + age_sc*man + 
                   educ_sc*man + inc_sc*man + pid7_dem*man, 
                 cem_dat_black, weights = mat_black$w)
summary(m_gend_blk)

# White respondents
m_gend_wht <- lm(cjs_att ~ pol_quality*man + pol.stop*man + rr_sc*man + peer.felony*man + age_sc*man + 
                   educ_sc*man + inc_sc*man + pid7_dem*man, 
                 cem_dat_white, weights = mat_white$w)
summary(m_gend_wht)


# Figure 7 --------------------------------------------------------------------
# * Race Figure Data Sets -------------------------------------------------
blk_samp_qual <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 1), na.rm = T)[[1]],
                            man = names(svytable(~man, subset(d.bw.stk, black == 1)))[which.max(svytable(~man, subset(d.bw.stk, black == 1)))],
                            educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 1)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 1)))],
                            inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 1)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 1)))],
                            pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 1)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 1)))],
                            pol.stop = 1/3,
                            pol_quality = c(0, .5, 1),
                            peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 1)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 1)))],
                            rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 1), na.rm = T)[[1]],
                            black = 1,
                            cjs_att = NA)
blk_samp_qual <- as.data.frame(apply(blk_samp_qual, 2, as.numeric))

blk_samp_qual$pol_quality <- factor(blk_samp_qual$pol_quality,
                                    labels = c("Positive", "Mixed", "Negative"))


wht_samp_qual <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 0), na.rm = T)[[1]],
                            man = names(svytable(~man, subset(d.bw.stk, black == 0)))[which.max(svytable(~man, subset(d.bw.stk, black == 0)))],
                            educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 0)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 0)))],
                            inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 0)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 0)))],
                            pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 0)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 0)))],
                            pol.stop = 1/3,
                            pol_quality = c(0, .5, 1),
                            peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 0)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 0)))],
                            rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 0), na.rm = T)[[1]],
                            black = 0,
                            cjs_att = NA)
wht_samp_qual <- as.data.frame(apply(wht_samp_qual, 2, as.numeric))
wht_samp_qual$pol_quality <- factor(wht_samp_qual$pol_quality,
                                    labels = c("Positive", "Mixed", "Negative"))


# * Race Figure Predictions ----------------------------------------------------------
# predicted values
blk_qual <- predict(m_race, blk_samp_qual, 
                    se.fit = TRUE)
wht_qual <- predict(m_race, wht_samp_qual, 
                    se.fit = TRUE)

# store in data frame
eval_preds <- data.frame(var = "pol_quality",
                         fit = c(blk_qual$fit, 
                                 wht_qual$fit),
                         se = c(blk_qual$se.fit, 
                                wht_qual$se.fit),
                         black = c(rep(1, 3), rep(0, 3)),
                         x = c(seq(0, 1, length.out = 3), seq(0, 1, length.out = 3)))

eval_preds$black <- factor(eval_preds$black,
                           labels = c("White", "Black"))

# * Race Facet -------------------------------------------------------------
eval_preds$x <- factor(eval_preds$x,
                       labels = c("Positive\n(N=1014)", "Mixed\n(N=394)", "Negative\n(N=320)"))
# Figure
plot_race <- ggplot(eval_preds, aes(x = x, y = fit, group = black)) +
  geom_ribbon(aes(ymin = fit - 1.39*se, ymax = fit + 1.39*se, fill = black),
              alpha = .6) +
  geom_line(size = 1, aes(color = black)) +
  scale_color_manual(values = c("grey", "black")) +
  scale_fill_manual(values = c("grey", "black")) +
  theme_bw() +
  scale_y_continuous(labels = ff_num(zero = 0, digits = 2)) +
  labs(x = "Contact Quality", y = "Predicted Carceral State Attitude") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")
plot_race


# * Gender Figure Data Sets -------------------------------------------------
blk_man_samp_qual <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 1 & man == 1), na.rm = T)[[1]],
                                man = 1,
                                educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 1 & man == 1)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 1 & man == 1)))],
                                inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 1 & man == 1)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 1 & man == 1)))],
                                pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 1 & man == 1)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 1 & man == 1)))],
                                pol.stop = 1/3,
                                pol_quality = c(0, .5, 1),
                                peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 1 & man == 1)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 1 & man == 1)))],
                                rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 1 & man == 1), na.rm = T)[[1]],
                                black = 1,
                                cjs_att = NA)
blk_man_samp_qual <- as.data.frame(apply(blk_man_samp_qual, 2, as.numeric))
blk_man_samp_qual$pol_quality <- factor(blk_man_samp_qual$pol_quality,
                                        labels = c("Positive", "Mixed", "Negative"))


blk_woman_samp_qual <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 1 & man == 0), na.rm = T)[[1]],
                                  man = 0,
                                  educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 1 & man == 0)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 1 & man == 0)))],
                                  inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 1 & man == 0)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 1 & man == 0)))],
                                  pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 1 & man == 0)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 1 & man == 0)))],
                                  pol.stop = 1/3,
                                  pol_quality = c(0, .5, 1),
                                  peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 1 & man == 0)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 1 & man == 0)))],
                                  rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 1 & man == 0), na.rm = T)[[1]],
                                  black = 1,
                                  cjs_att = NA)
blk_woman_samp_qual <- as.data.frame(apply(blk_woman_samp_qual, 2, as.numeric))
blk_woman_samp_qual$pol_quality <- factor(blk_woman_samp_qual$pol_quality,
                                          labels = c("Positive", "Mixed", "Negative"))


wht_man_samp_qual <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 0 & man == 1), na.rm = T)[[1]],
                                man = 1,
                                educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 0 & man == 1)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 0 & man == 1)))],
                                inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 0 & man == 1)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 0 & man == 1)))],
                                pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 0 & man == 1)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 0 & man == 1)))],
                                pol.stop = 1/3,
                                pol_quality = c(0, .5, 1),
                                peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 0 & man == 1)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 0 & man == 1)))],
                                rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 0 & man == 1), na.rm = T)[[1]],
                                black = 0,
                                cjs_att = NA)
wht_man_samp_qual <- as.data.frame(apply(wht_man_samp_qual, 2, as.numeric))
wht_man_samp_qual$pol_quality <- factor(wht_man_samp_qual$pol_quality,
                                        labels = c("Positive", "Mixed", "Negative"))


wht_woman_samp_qual <- data.frame(age_sc = svymean(~age_sc, subset(d.bw.stk, black == 0 & man == 0), na.rm = T)[[1]],
                                  man = 0,
                                  educ_sc = names(svytable(~educ_sc, subset(d.bw.stk, black == 0 & man == 0)))[which.max(svytable(~educ_sc, subset(d.bw.stk, black == 0 & man == 0)))],
                                  inc_sc = names(svytable(~inc_sc, subset(d.bw.stk, black == 0 & man == 0)))[which.max(svytable(~inc_sc, subset(d.bw.stk, black == 0 & man == 0)))],
                                  pid7_dem = names(svytable(~pid7_dem, subset(d.bw.stk, black == 0 & man == 0)))[which.max(svytable(~pid7_dem, subset(d.bw.stk, black == 0 & man == 0)))],
                                  pol.stop = 1/3,
                                  pol_quality = c(0, .5, 1),
                                  peer.felony = names(svytable(~peer.felony, subset(d.bw.stk, black == 0 & man == 0)))[which.max(svytable(~peer.felony, subset(d.bw.stk, black == 0 & man == 0)))],
                                  rr_sc = svymean(~rr_sc, subset(d.bw.stk, black == 0 & man == 0), na.rm = T)[[1]],
                                  black = 0,
                                  cjs_att = NA)
wht_woman_samp_qual <- as.data.frame(apply(wht_woman_samp_qual, 2, as.numeric))
wht_woman_samp_qual$pol_quality <- factor(wht_woman_samp_qual$pol_quality,
                                          labels = c("Positive", "Mixed", "Negative"))


# * Gender Figure Predictions ----------------------------------------------------------
# predicted values
wht_man_direct <- predict(m_gend_wht, wht_man_samp_qual, 
                          se.fit = TRUE)
wht_wom_direct <- predict(m_gend_wht, wht_woman_samp_qual, 
                          se.fit = TRUE)

# store in data frame
pred_direct<- data.frame(var = "direct_exp",
                         fit = c(wht_man_direct$fit, 
                                 wht_wom_direct$fit),
                         se = c(wht_man_direct$se.fit, 
                                wht_wom_direct$se.fit),
                         man = c(rep(1, 3), rep(0, 3)),
                         x = c(seq(0, 1, length.out = 3), seq(0, 1, length.out = 3)))
eval_preds_wht <- rbind(pred_direct)
eval_preds_wht$man <- factor(eval_preds_wht$man,
                             labels = c("Woman", "Man"))


# predicted values
blk_man_direct <- predict(m_gend_blk, blk_man_samp_qual, 
                          se.fit = TRUE)
blk_wom_direct <- predict(m_gend_blk, blk_woman_samp_qual, 
                          se.fit = TRUE)

# store in data frame
pred_direct<- data.frame(var = "direct_exp",
                         fit = c(blk_man_direct$fit, 
                                 blk_wom_direct$fit),
                         se = c(blk_man_direct$se.fit, 
                                blk_wom_direct$se.fit),
                         man = c(rep(1, 3), rep(0, 3)),
                         x = c(seq(0, 1, length.out = 3), seq(0, 1, length.out = 3)))
eval_preds_blk <- rbind(pred_direct)
eval_preds_blk$man <- factor(eval_preds_blk$man,
                             labels = c("Woman", "Man"))
eval_preds_wht$race <- "White"
eval_preds_blk$race <- "Black"

# combine
preds_comb <- rbind(eval_preds_wht, eval_preds_blk)


# * Gender Facet -------------------------------------------------------------
preds_comb$x2 <- NA
preds_comb$x2[which(preds_comb$race == "White")] <- c("Positive\n(N=755)", "Mixed\n(N=172)", "Negative\n(N=125)")
preds_comb$x2[which(preds_comb$race == "Black")] <- c("Positive\n(N=256)", "Mixed\n(N=165)", "Negative\n(N=160)")
preds_comb$x2 <- factor(preds_comb$x2,
                        levels = c("Positive\n(N=755)", "Positive\n(N=256)", "Mixed\n(N=172)",
                                   "Mixed\n(N=165)", "Negative\n(N=125)", "Negative\n(N=160)"))
preds_comb$man <- ifelse(preds_comb$man == "Man", "Men", "Women")

# Figure
plot_gend <- ggplot(preds_comb, aes(x = x2, y = fit, group = man)) +
  geom_ribbon(aes(ymin = fit - 1.39*se, ymax = fit + 1.39*se, fill = man),
              alpha = .6) +
  geom_line(size = 1, aes(color = man)) +
  scale_color_manual(values = c("black", "grey")) +
  scale_fill_manual(values = c("black", "grey")) +
  scale_y_continuous(labels = ff_num(zero = 0, digits = 2)) +
  facet_grid(~race, scales = "free_x") +
  theme_bw() +
  labs(x = "Contact Quality", y = "Predicted Carceral State Attitude") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size = 14),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.text = element_text(size = 14),
        legend.position = "bottom")
plot_gend

# * Combine and Plot -------------------------------------------------------------

(plot_race / plot_gend) + plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 18))
# ggsave("figure7.pdf", width = 8, height = 8)

